###########################################################################
## Challenging Encounters and Within-Physician Practice Variability
## Script 4- Analysis
## Last update: 29-1-23
###########################################################################

## Setup
###########################################################################

# Packages
library(dplyr)   
library(xtable)   
library(readxl)   
library(lmtest)    
library(lubridate)
library(miceadds)
library(plm)
library(stargazer)

# Directory
wd <- "."
setwd(wd)

start_file <- proc.time() #Keep File Starting Time
###########################################################################

# Functions
###########################################################################

# Adjust coeftest function to include clustered SEs (c is a plm regression object)
coeftest_cluster <- function(x){
    coeftest(x, vcov. = vcovHC(x, type = "HC0", cluster = "group"))
  }

# Dynamic diff-in-diff
dynamic_dd <- function(outcome, min1, max1, Label) {
  # Arguments -----------------------
  # Outcome: outcome variable
  # min1: lower limit of y axis
  # max1: higher limit of y axis
  # Label:  label for y
  # ---------------------------------
  
  # Estimate dynamic DD models with the three different comparison groups
  Reg1 <-plm(formula = paste(outcome,
                             "~TBI_Pre_6+TBI_Pre_5+TBI_Pre_4+TBI_Pre_3+TBI_Pre_2+TBI_1+TBI_2+TBI_3+TBI_4+TBI_5+TBI_6+TBI_7+TBI_8+TBI_9+BI_Pre_6+BI_Pre_5+BI_Pre_4+BI_Pre_3+BI_Pre_2+BI_1+BI_2+BI_3+BI_4+BI_5+BI_6+BI_7+BI_8+BI_9+Treat",
                             sep = ""),data = Vis[Vis$Treat == 1 | Vis$Control==1 ,],index = c("Event_ID"),model = "within")
  
  Reg2 <-plm(formula = paste(outcome,
                             "~TBI_Pre_6+TBI_Pre_5+TBI_Pre_4+TBI_Pre_3+TBI_Pre_2+TBI_1+TBI_2+TBI_3+TBI_4+TBI_5+TBI_6+TBI_7+TBI_8+TBI_9+BI_Pre_6+BI_Pre_5+BI_Pre_4+BI_Pre_3+BI_Pre_2+BI_1+BI_2+BI_3+BI_4+BI_5+BI_6+BI_7+BI_8+BI_9+Treat",
                             sep = ""),data = Vis[Vis$Treat == 1 | Vis$Control==2 ,],index = c("Event_ID"),model = "within")
  
  Reg3 <-plm(formula = paste(outcome,
                             "~TBI_Pre_6+TBI_Pre_5+TBI_Pre_4+TBI_Pre_3+TBI_Pre_2+TBI_1+TBI_2+TBI_3+TBI_4+TBI_5+TBI_6+TBI_7+TBI_8+TBI_9+BI_Pre_6+BI_Pre_5+BI_Pre_4+BI_Pre_3+BI_Pre_2+BI_1+BI_2+BI_3+BI_4+BI_5+BI_6+BI_7+BI_8+BI_9+Treat",
                             sep = ""),data = Vis[Vis$Treat == 1 | Vis$Control==3 ,],index = c("Event_ID"),model = "within")
  
  # Keep estimates
  coefcl1 <- coeftest_cluster(Reg1)
  coefcl2 <- coeftest_cluster(Reg2)
  coefcl3 <- coeftest_cluster(Reg3)
  
  coef1 <- c(coefcl1[1:5, 1], 0, NA, coefcl1[6:14, 1])
  se1 <- c(coefcl1[1:5, 2], 0, NA, coefcl1[6:14, 2])
  
  coef2 <- c(coefcl2[1:5, 1], 0, NA, coefcl2[6:14, 1])
  se2 <- c(coefcl2[1:5, 2], 0, NA, coefcl2[6:14, 2])
  
  coef3 <- c(coefcl3[1:5, 1], 0, NA, coefcl3[6:14, 1])
  se3 <- c(coefcl3[1:5, 2], 0, NA, coefcl3[6:14, 2])
  
  
  R1 <- coef1
  R1_Low <- coef1 - 1.645 * se1
  R1_Hig <- coef1 + 1.645 * se1
  
  R2 <- coef2
  R2_Low <- coef2 - 1.645 * se2
  R2_Hig <- coef2 + 1.645 * se2
  
  R3 <- coef3
  R3_Low <- coef3 - 1.645 * se3
  R3_Hig <- coef3 + 1.645 * se3
  
  # Plot DD estimates
  pdf(file.path(pd,"outputs","Figure_1b.pdf"))
  plot(R1 ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1),cex = 0.75,col = scales::alpha("blue1",0.3),
       bty = "n",xlim = c(1, 16),xlab = "Visits Relative to Index Event",
       ylab = Label, pch = 15,type = "p",main = "")
  par(new = T)
  
  plot(R3 ~ c(1:16+0.1),xaxt = "n",yaxt = "n",ylim = c(min1, max1),cex = 0.75,col = "blue2",
       bty = "n",xlim = c(1, 16),xlab = "",
       ylab ="",pch = 16,type = "l",main = "")
  par(new = T)
  plot(R3 ~ c(1:16+0.1),xaxt = "n",yaxt = "n",ylim = c(min1, max1),cex = 0.75,col = "blue2",
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 16,type = "p",main = "")
  par(new = T)
  
  plot(R2 ~ c(1:16+0.2),xaxt = "n",yaxt = "n",ylim = c(min1, max1),cex = 0.75,col = scales::alpha("blue3",0.3),
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 17,type = "p",main = "")
  
  axis(1,at = seq(1, 16, by = 1),labels = seq(-6, 9, by = 1)*2,cex.axis = 0.75,las = 1,pos = min1)
  axis(2,at = seq(min1, max1 + 1, by = (max1 - min1) / 8),cex.axis = 0.75,las = 2,pos = 0.5)
  abline(h = max1, col = scales::alpha("darkgrey", 0.25))
  abline(h = 0.5 * max1, col = scales::alpha("darkgrey", 0.25))
  abline(h = 0.5 * min1, col = scales::alpha("darkgrey", 0.25))
  abline(h = 0, col = scales::alpha("black", 0.75))
  lines(x = c(7, 7),y = c(min1, 100),col = "gray8",lty = 5)
  abline(h = min1,col = scales::alpha("black"),lwd = 1.1)
  
  arrows(1:16,R1_Low,1:16,R1_Hig,lty = 1,lwd = 1,col = scales::alpha("blue1", 0.3),angle = 90,
         length = 0.03,code = 3)
  
  arrows(1:16+0.1,R3_Low,1:16+0.1,R3_Hig,lty = 1,lwd = 1,col = scales::alpha("blue2", 0.8),angle = 90,
         length = 0.03,code = 3)
  
  arrows(1:16+0.2,R2_Low,1:16+0.2,R2_Hig,lty = 1,lwd = 1,col = scales::alpha("blue3", 0.3),angle = 90,
         length = 0.03,code = 3)
  
  text(x =8.5, y = max1-(max1-min1)/20,labels = "Index Event", col = "gray8")
  
  dev.off()
  
  
  # Flexible time-trends
  
  coef1 <- c(coefcl1[1:5, 1], 0, NA, coefcl1[6:14, 1])
  se1 <- c(coefcl1[1:5, 2], 0, NA, coefcl1[6:14, 2])
  
  coef2 <- c(coefcl2[1:5, 1], 0, NA, coefcl2[6:14, 1])
  se2 <- c(coefcl2[1:5, 2], 0, NA, coefcl2[6:14, 2])
  
  coef3 <- c(coefcl3[1:5, 1], 0, NA, coefcl3[6:14, 1])
  coef3C <- c(coefcl3[15:19, 1], 0, NA, coefcl3[20:28, 1])
  
  
  R1 <- coef1
  
  R2 <- coef2
  
  R3 <- coef3
  
  
  R1T <- coef3C + coef3
  R1C1 <- coef3C + coef3 - coef1
  R1C2 <- coef3C + coef3 - coef2
  R1C3 <- coef3C
  
  pdf(file.path(pd,"outputs","Figure_1a.pdf"))
  plot(R1T ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = "darkred",
       bty = "n",xlim = c(1, 16),xlab = "Visits Relative to Index Event",ylab = Label,pch = 16,
       main = "",type = "l")
  
  par(new=T)
  plot(R1T ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = "darkred",
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 16,
       main = "",type = "p")
  
  axis(1,at = seq(1, 16, by = 1),labels = seq(-6, 9, by = 1)*2,cex.axis = 0.75,las = 1,pos = min1)
  axis(2,at = seq(min1, max1 + 0.5, by = (max1 - min1) / 8),cex.axis = 0.75,las = 2,pos = 0.8)
  
  abline(h = seq(min1,max1,(max1 - min1) / 4),
         col = scales::alpha("darkgrey", 0.25))
  abline(h = min1,
         col = scales::alpha("black"),
         lwd = 1.1)
  lines(
    x = c(7, 7),
    y = c(min1, 100),
    col = "gray8",
    lty = 5
  )
  par(new = TRUE)
  plot(R1C1 ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = scales::alpha("blue1",0.3),
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 18,lty = 4,main = "",type = "l")
  
  par(new = TRUE)
  plot(R1C2 ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = scales::alpha("blue2",0.3),
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 17,lty = 3,main = "",type = "l")
  
  par(new = TRUE)
  plot(R1C3 ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = "blue3",
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 15,lty = 2,main = "",type = "l")
  
  par(new = TRUE)
  plot(R1C1 ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = scales::alpha("blue1",0.3),
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 18,lty = 4,main = "",type = "p")
  
  par(new = TRUE)
  plot(R1C2 ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = scales::alpha("blue2",0.3),
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 17,lty = 3,main = "",type = "p")
  
  par(new = TRUE)
  plot(R1C3 ~ c(1:16),xaxt = "n",yaxt = "n",ylim = c(min1, max1 + 0.001),cex = 0.75,col = "blue3",
       bty = "n",xlim = c(1, 16),xlab = "",ylab = "",pch = 15,lty = 2,main = "",type = "p")
  
  
  lines(x = c(6,8),y = c(R1T[6],R1T[8]),col = "darkred", lty = 1)
  lines(x = c(6,8),y = c(R1C1[6],R1C1[8]),col = scales::alpha("blue1",0.3), lty = 2)
  lines(x = c(6,8),y = c(R1C2[6],R1C2[8]),col = scales::alpha("blue2",0.3), lty = 3)
  lines(x = c(6,8),y = c(R1C3[6],R1C3[8]),col = "blue3", lty = 4)
  
  text(x =11.9, y = 1.7,labels = "Difficult", col = "darkred",cex = 1.5)
  text(x =12.1, y = 1.25,labels = "Index Case", col = "darkred",cex = 1.5)
  text(x =9.2, y = -1.2,labels = "Comparison", col = "blue3",cex = 1.5)
  text(x =9.3, y = -1.65,labels = "Index Cases", col = "blue3",cex = 1.5)
  
  dev.off()
  
}

# diff-in-diff
DD_main <- function(outcomes, names, CON=3, v1=8 ,name){
  # Arguments -----------------------
  # Outcomes: outcome variables
  # names: outcome names
  # CON: comparison group (3 is the main and the default)
  # v1: sample restriction on the number of visits before/after the index case (default is 8)
  # name: Table's name
  # ---------------------------------
  
  # Restrict the sample (comparison group and window of v1 visits before/after)
  V <- Vis
  Vis <- V[abs(V$Visits_From_Event) < v1+1 & (V$Treat==1 | V$Control == CON),]
  
  
  # Loop over outcomes
  Table <- NULL
  i <- 0
  for(outcome in outcomes){
    i <- i+1
    
    # Estimate a DD model
    if(outcome=="visit_length"){ # In the case of visit duration, we exclude the visit before the index case due to a measurement error
      Post2 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat", sep = ""),
                   data = Vis[Vis$Visits_From_Event!=(-1),],index = c("Real_Event_ID"),model = "within")
    }else{
      Post2 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat", sep = ""),
                   data = Vis,index = c("Real_Event_ID"),model = "within")
    }
    
    # Keep estimates
    coefcl2 <- coeftest_cluster(Post2)
    est2 <- round(coefcl2[1, 1], 2)
    sd2 <- round(coefcl2[1, 2], 2)
    t2 <- round(coefcl2[1, 3], 2)
    N2 <- length(Post2$residuals)
    CL1 <- length(unique(Vis$Real_Event_ID))
    Baseline <- mean(V[[outcome]][V$Visits_From_Event < 0 & (V$Treat==1 | V$Control == CON)])
    
    # Export the specific regression (with full details)
    write.csv(rbind(coefcl2,N2,CL1),file.path(pd,"outputs","other",
                                              paste("DD_all_coefs", CON, outcome,"v1",v1, ".csv", sep = "_")))
    
    # Update table
    Table <- rbind(Table,as.data.frame(list("Outcome" = c(names[i],"",""),
                                            "Baseline" = c(Baseline,NA,NA),
                                            "Estimate" = c(est2,sd2,NA))))
    
  }
  
  # Export the full table (paper's version)
  table <- print(xtable(Table))
  write.table(table,file.path(pd,"outputs",paste(name,".txt", sep = "_")))

  # Restore
  Vis <- V
}

# diff-in-diff with control variables
DD_main_controls <- function(outcomes, names, CON=3, v1=8, name) {
  # Arguments -----------------------
  # Outcomes: outcome variables
  # names: outcome names
  # CON: comparison group (3 is the main and the default)
  # v1: sample restriction on the number of visits before/after the index case (default is 8)
  # name: Table's name
  # ---------------------------------
  
  # Define the controls
  C2 <- "+TIA+CVD+Diabetic+Hashmana+Cancer_Old+age+age_sq+gender+Hospitalizations_Last_6_Months+Days_Hospitalizations_Last_6_Months+Visits_Last_6_Months+Visits_Last_6_Months_PCP"
  
  # Restrict the sample (comparison group and window of v1 visits before/after)
  V <- Vis
  Vis <- V[abs(V$Visits_From_Event) < v1+1 & (V$Treat==1 | V$Control == CON),]
  
  # Loop over outcomes
  Table <- NULL
  i <- 0
  for(outcome in outcomes){
    i <- i+1
    
    # Estimate a DD model
    if(outcome=="visit_length"){ # In the case of visit duration, we exclude the visit before the index case due to a measurement error
      Post2 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat",C2, sep = ""),
                   data = Vis[Vis$Visits_From_Event!=(-1),],index = c("Real_Event_ID"),model = "within")
    }else{
      Post2 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat",C2, sep = ""),
                   data = Vis,index = c("Real_Event_ID"),model = "within")
    }
    
    # Keep estimates
    coefcl2 <- coeftest_cluster(Post2)
    est2 <- round(coefcl2[1, 1], 2)
    sd2 <- round(coefcl2[1, 2], 2)
    t2 <- round(coefcl2[1, 3], 2)
    N2 <- length(Post2$residuals)
    CL1 <- length(unique(Vis$Real_Event_ID))
    Baseline <- mean(V[[outcome]][V$Visits_From_Event < 0 & (V$Treat==1 | V$Control == CON)])
    
    # Export the specific regression (with full details)
    write.csv(rbind(coefcl2,N2,CL1),file.path(pd,"outputs","other",
                                              paste("DD_controls_all_coefs", CON, "v1",v1,".csv", sep = "_")))
    
    # Update table
    Table <- rbind(Table,as.data.frame(list("Outcome" = c(names[i],"",""),
                                            "Baseline" = c(Baseline,NA,NA),
                                            "Estimate" = c(est2,sd2,NA))))
    
  }
  
  # Export the full table (paper's version)
  table <- print(xtable(Table))
  write.table(table,file.path(pd,"outputs",paste(name,".txt", sep = "_")))
  
  
  Vis <- V
}

# diff-in-diff, compare the two alternative comparison groups
DD_main_alternative <- function(outcomes, names,v1=8,name){
  # Arguments -----------------------
  # Outcomes: outcome variables
  # names: outcome names
  # v1: sample restriction on the number of visits before/after the index case (default is 8)
  # name: Table's name
  # ---------------------------------
  
  # Restrict the sample (window of v1 visits before/after)
  V <- Vis
  Vis <- V[abs(V$Visits_From_Event) < v1+1,]
  
  
  # Loop over outcomes
  Table <- NULL
    i <- 0
  for(outcome in outcomes){
    i <- i+1
    
    # Estimate
    if(outcome=="visit_length"){# In the case of visit duration, we exclude the visit before the index case due to a measurement error
      Post1 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat", sep = ""),
                   data = Vis[Vis$Visits_From_Event!=(-1) & Vis$Control %in% c(0,1),],index = c("Real_Event_ID"),model = "within")
      Post2 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat", sep = ""),
                   data = Vis[Vis$Visits_From_Event!=(-1) & Vis$Control %in% c(0,2),],index = c("Real_Event_ID"),model = "within")
    }else{
      Post1 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat", sep = ""),
                   data = Vis[Vis$Control %in% c(0,1),],index = c("Real_Event_ID"),model = "within")
      Post2 <- plm(formula = paste(outcome, "~Treat_Post+Post+Treat", sep = ""),
                   data = Vis[Vis$Control %in% c(0,2),],index = c("Real_Event_ID"),model = "within")
    }
    
    # Keep estimates 
    coefcl1 <- coeftest_cluster(Post1)
    est1 <- round(coefcl1[1, 1], 2)
    sd1 <- round(coefcl1[1, 2], 2)
    N1 <- length(Post1$residuals)
    CL1 <- length(unique(Vis$Real_Event_ID[Vis$Control == 1]))
    
    coefcl2 <- coeftest_cluster(Post2)
    est2 <- round(coefcl2[1, 1], 2)
    sd2 <- round(coefcl2[1, 2], 2)
    N2 <- length(Post2$residuals)
    CL2 <- length(unique(Vis$Real_Event_ID[Vis$Control == 2]))
    
    Baseline <- mean(V[[outcome]][V$Visits_From_Event < 0 & (V$Treat==1 | V$Control == 3)])
    
    # Export the specific regression (with full details)
    write.csv(rbind(coefcl2,N2,CL2),file.path(pd,"outputs","other",
                                              paste("DD_all_coefs_alternative", 2, outcome,"v1",v1, ".csv", sep = "_")))
    write.csv(rbind(coefcl1,N1,CL1),file.path(pd,"outputs","other",
                                              paste("DD_all_coefs_alternative", 1, outcome,"v1",v1, ".csv", sep = "_")))
    
    # Update table
    Table <- rbind(Table,as.data.frame(list("Outcome" = c(names[i],"",""),
                                            "Baseline" = c(Baseline,NA,NA),
                                            "Estimate (1)" = c(est1,sd1,NA),
                                            "Estimate (2)" = c(est2,sd2,NA))))
    
  }
  
  # Export the paper's version
  table <- print(xtable(Table))
  write.table(table,file.path(pd,"outputs",paste(name,".txt", sep = "_")))
  
  # Restore
  Vis <- V
}

# triple-diffs
DDD <- function(outcome,Z, names, CON=3,v1=8,name){
  # Arguments -----------------------
  # Outcomes: outcome variables
  # Z: heterogeneity dimensions
  # names: outcome names
  # CON: comparison group (3 is the main and the default)
  # v1: sample restriction on the number of visits before/after the index case (default is 8)
  # name: Table's name
  # ---------------------------------
  
  # Restrict the sample (comparison group and window of v1 visits before/after)
  V <- Vis
  Vis <- V[abs(V$Visits_From_Event) < v1+1 & (V$Treat==1 | V$Control == CON),]
  
  # Loop over outcomes
  Table <- NULL
  i <- 0
  for(z in Z){
    i <- i+1
    
    # Estimate DDD model
    Vis$Char <- Vis[[z]]
    Vis$Post_Char <- Vis$Char * Vis$Post
    Vis$Treat_Char <- Vis$Char * Vis$Treat
    Vis$Treat_Post_Char <- Vis$Char * Vis$Post * Vis$Treat
    Post2 <-plm(formula = paste(outcome,"~Treat_Post_Char+Treat_Post+Treat_Char+Post_Char+Treat+Post+Char",sep = ""),
                data = Vis,index = c("Real_Event_ID"),model = "within")
    
    # Keep estimates
    coefcl2 <- coeftest_cluster(Post2)
    coef2 <- coeftest_cluster(Post2)
    est2 <- round(coef2[1, 1], 4)
    sd2 <- round(coef2[1, 2], 4)
    t2 <- round(coef2[1, 3], 4)
    N2 <- length(Post2$residuals)
    CL1 <- length(unique(Vis$Real_Event_ID))
    
    # Export estimates
    write.csv(rbind(coefcl2,N2,CL1),file.path(pd,"outputs","other",
                                              paste("DDD_all_coefs", "v1",v1,CON, outcome,"_",z, ".csv", sep = "_")))
    
    # Export the implied effects (for each value of z)
    X <- unique(Vis[[z]])
    X <- X[order(X)]
    Effect <- coef2[2,1]+coef2[1,1] * X
    VC <- vcovHC(Post2, type = "HC0", cluster = "group")
    SE <- sqrt(VC[2,2] + VC[1,1]*X^2 + 2*X*VC[1,2])
    write.csv(rbind(Effect,SE),
              file =file.path(pd,"outputs","other",paste("DDD_effects_",outcome,"_",z,"v1",v1,".csv",sep="")))
    
    
    # Update table
    Table <- rbind(Table,as.data.frame(list("Outcome" = c(names[i],"",""),
                                            "No" = c(Effect[1],SE[1],NA),
                                            "Yes" = c(Effect[2],SE[2],NA),
                                            "Difference" = c(est2,sd2,NA))))
    
  }
  
  # Export the full table (paper's version)
  table <- print(xtable(Table))
  write.table(table,file.path(pd,"outputs",paste(name,".txt", sep = "_")))
  
  # Restore
  Vis <- V
}

# triple-diffs with the estimated propensity score
PS_DDD <- function(outcome,CON=3,v1=8){
  # Arguments -----------------------
  # Outcome: the estimated propensity scores
  # CON: comparison group (3 is the main and the default)
  # v1: sample restriction on the number of visits before/after the index case (default is 8)
  # ---------------------------------
  
  # Preserve
  V <- Vis
  
  # Mark the propensity score
  Vis$px <-Vis[[outcome]]
  
  # Normalize 
  Vis$R <-Vis[[outcome]]
  c <- quantile(Vis$px ,probs=seq(0.01,1,0.01))
  Vis$R <- 1
  for(i in 100:1){
    Vis$R[Vis[[outcome]]<c[i]] <- (i-1)*0.01
  }
  Vis$R[Vis[[outcome]]==c[i]] <- 0
  
  # Calibration plot's (Figure A5) data
  plot <- summarise(group_by(Vis,R,Post),mean_px=mean(px),mean_diag=mean(Diagnostic_Test),total = n())
  write.csv(plot,"Emotions/outputs/Calibration_plot.csv")

  # Plot
  pdf("Emotions/outputs/Figure_A5.pdf")
  plot(plot$mean_diag[plot$Post==0] ~ plot$mean_px[plot$Post==0], bty = "n", xlab = "Average Prediction", xlim = c(0,0.6), ylim = c(0,0.6),
       ylab = "Average Rate of Testing", xaxt = "n",yaxt = "n", pch = 1, col = "blue3", cex = 1)
  abline(0,1, col = "darkred")
  axis(1,seq(0,1,0.1),pos=0)
  axis(2,seq(0,1,0.1),pos=0)
  dev.off()
  
  # Estimate a DDD Model
  Vis$R_Treat_Post <-Vis$R*Vis$Treat_Post
  Vis$R_Post <-Vis$R*Vis$Post
  Vis$R_Treat <-Vis$R*Vis$Treat
  Vis <- Vis[abs(Vis$Visits_From_Event) < v1+1 & (Vis$Treat==1 | Vis$Control == CON),]
  DDD <- plm(formula=paste("Diagnostic_Test ~ R+ R_Treat_Post+Treat_Post+Treat+Post+R_Treat+R_Post",sep=""),data= Vis,index = c("Real_Event_ID"), model = "within")

  # Keep estimates
  coef <- coeftest_cluster(DDD)

  # Calculate the implied effects
  X <- seq(0,1,0.01)
  Effect <- coef[3,1]+coef[2,1] * X
  CL <- vcovHC(DDD, type = "HC0", cluster = "group")
  SE <- sqrt(CL[3,3] + CL[2,2] * X^2 + 2*X*CL[3,2])

# Plot Figure 2
Effect <- as.numeric(Effect)
Min <- Effect - 1.645 * as.numeric(SE)
Max <- Effect + 1.645 * as.numeric(SE)

pdf(file.path(pd,"outputs","Figure_2.pdf"))
plot(Effect ~ c(seq(0,1,0.01)), type = "l", bty = "n", xaxt = "n", ylab = "", 
     xlab = "", col ="blue3", yaxt = "n", xlim = c(0.02,1), 
     ylim = c( -4, 4.2))

mtext("Propensity Score for Testing", side=1, line=2.2, cex=1.5)
mtext("Estimated Effect on Rate of Testing (p.p.)", side=2, line=2.5, cex=1.3)

abline(h=0)
abline(h=seq(-2,4,2), col = "grey88")
par(new=T)
plot(Min ~ c(seq(0,1,0.01)), type = "l", bty = "n", xaxt = "n", ylab = "", 
     xlab = "",yaxt = "n",lty = 2, col ="darkred", xlim = c(0.02,1),
     ylim = c( -4, 4.2))
par(new=T)
plot(Max ~ c(seq(0,1,0.01)), type = "l", bty = "n", xaxt = "n", ylab = "", 
     xlab = "",yaxt = "n", lty = 2, col ="darkred", xlim = c(0.02,1),
     ylim = c( -4, 4.2))
axis(1,seq(-1,1,0.1), pos = -4,cex.axis=1.3)
axis(2,seq(-4,5,1), pos = 0.0,cex.axis=1.3)
dev.off()

  
}

###########################################################################

## Import data and sample restrictions
###########################################################################

# Import data:
Vis <-readRDS(file.path(pd,"data","data_final_with_prediction.R"))

# Add analysis columns
Vis$Post <- as.numeric(Vis$Visits_From_Event > 0)
Vis$Treat_Post <- as.numeric(Vis$Treat * Vis$Post)

# Keep cases with at least 1 comparison case
sum <- summarise(
  group_by(Vis, Real_Event_ID),
  Events = sum(Event),
  Treat = sum(Event_Treat),
  Control0 = sum(Event * Control==0),
  Control1 = sum(Event * Control==1),
  Control2 = sum(Event * Control==2),
  Control3 = sum(Event * Control==3))
sum <- sum[sum$Treat > 0 & (sum$Events > 1), ]
sum$w <- 1
Vis <- merge(Vis, sum[, c("Real_Event_ID", "w")], by = "Real_Event_ID")

# Sample of 12 Visits before and 18 after the index case
Vis <- Vis[Vis$Visits_From_Event > -13 & Vis$Visits_From_Event < 19, ]

###########################################################################

# Appendix figures A1-A4
###########################################################################

# Index cases
Events <- Vis[Vis$Event == 1 , ]
Events <- Events[order(Events$doctor_id, Events$date),]

# Exposure to difficult cases (Figure A1)
Do <- summarise(group_by(Events, doctor_id), total = sum(Treat))
Doc <- summarise(group_by(Do, total), Docs = n())
Doc <- Doc[Doc$total != 0, ]
pdf(file.path(pd,"outputs","Figure_A1.pdf"))
barplot(
  Doc$Docs,
  ylab = "Number of Physicians",xaxt="n",yaxt="n",
  xlab = "Number of Difficult Cases",col="darkred",
  names.arg = 1:length(t(Doc[, 1]))
)
axis(1,at=seq(0.7,0.7+2.4*49,2.4),labels=seq(1,100,2),pos=0, cex.axis = 0.8)
axis(2,seq(0,140,20), cex.axis = 0.8)
abline(h=0)
dev.off()

# Add indicator
Vis <- merge(Vis,Do,by="doctor_id")
Vis$Exposure_5<-as.numeric(Vis$total>median(Do$total))

# Cases by their hour (Figure A2)
Events$Time_Fixed <- floor(Events$Time_Fixed * 2) / 2
Time <-summarise(group_by(Events, Time_Fixed),All_Events = n(),
                 Real_Events = sum(Treat),
                 Comp_Events1 = sum(Control==1),
                 Comp_Events2 = sum(Control==2),
                 Comp_Events3 = sum(Control==3))

Time_All <- readRDS(file.path(pd,"data","Visits_by_hour")) # Import all visits data

pdf(file.path(pd,"outputs","Figure_A2.pdf"))
plot(Time_All$Visits *100/ sum(Time_All$Visits) ~ Time_All$`floor(Visits$Time_Fixed * 2)/2`,
     ylim = c(0,10),type = "l",col = scales::alpha("green4", 0.6),pch = 8,xlab = "",ylab = "",lty=5,lwd=1.2,
     xlim = c(0, 24),xaxt = "n",yaxt = "n", bty = "n")
par(new = T)
plot(Time$Real_Events  *100/ sum(Time$Real_Events) ~ Time$Time_Fixed,lwd=1.2,
     ylim = c(0,10),lty = 1,type = "l",col = "darkred",pch = 16,xlab = "Time of Day",
     ylab = "Share of Cases (%)",xaxt = "n",xlim = c(0, 24), bty = "n")

par(new = T)
plot(Time$Comp_Events3 *100 / sum(Time$Comp_Events3) ~ Time$Time_Fixed,lwd=1.2,
     ylim = c(0,10),lty = 2,type = "l",col = "blue3",pch = 15,xlab = "",
     ylab = "",xlim = c(0, 24),xaxt = "n",yaxt = "n", bty = "n")
par(new = T)

plot(Time$Comp_Events2 *100 / sum(Time$Comp_Events2) ~ Time$Time_Fixed,lwd=1.2,
     ylim = c(0,10),lty = 3,type = "l",col = scales::alpha("blue2",0.2),pch = 17,
     xlab = "",ylab = "",xlim = c(0, 24),xaxt = "n",
     yaxt = "n", bty = "n")
par(new = T)

plot(Time$Comp_Events1 *100 / sum(Time$Comp_Events1) ~ Time$Time_Fixed,lwd=1.2,
     ylim = c(0,10),lty = 4,type = "l",col =  scales::alpha("blue1",0.2),pch = 18,
     xlab = "",ylab = "",xlim = c(0, 24),xaxt = "n",yaxt = "n", bty = "n")

legend(0,10,legend = c("Difficult Cases","Comparison Cases (Main)","Comparison Cases (Alternative I)","Comparison Cases (Alternative II)","All Cases"),
       lty = c(1,4,3,2,5), bty = "n", col = c("darkred","blue3","blue1","blue2","palegreen3"), cex = 0.75)
axis(1,seq(-4,28,4), pos =0)
dev.off()


# Cases by their month (Figure A3)
Monthly <-
  summarise(
    group_by(Events, year, month),
    All_Events = n(),
    Real_Events = sum(Treat),
    Comp_Events1 = sum(Control==1),
    Comp_Events2 = sum(Control==2),
    Comp_Events3 = sum(Control==3)
  )
Monthly_All <- readRDS(file.path(pd,"data","Visits_by_month"))

pdf(file.path(pd,"outputs","Figure_A3.pdf"))
plot(Monthly_All$Visits  *100/ sum(Monthly_All$Visits) ~ c(1:42),
     ylim = c(0, 5),type = "l",col = scales::alpha("green4", 0.6),pch = 8,
     xlab = "",xaxt = "n",ylab = "",lty = 5, bty = "n")
axis(1,seq(-5,48,6),
     labels=c("","Jul-12","Jan-13","Jul-13","Jan-14","Jul-14",
              "Jan-15","Jul-15","Jan-16"), pos =0 , cex.axis=0.85)
par(new = T)
plot(Monthly$Real_Events *100/ sum(Monthly$Real_Events) ~ c(1:42),
     ylim = c(0, 5),xaxt = "n",yaxt = "n",type = "l",col = "darkred",
     pch = 16,xlab = "",ylab = "Share of Cases (%)",
     lty = 1, bty = "n")

par(new = T)
plot(Monthly$Comp_Events3 *100 / sum(Monthly$Comp_Events3) ~ c(1:42),
     ylim = c(0, 5),type = "l",col = "blue3",pch = 15,xlab = "",
     ylab = "",lty = 2,xaxt = "n",yaxt = "n", bty = "n")

par(new = T)
plot(Monthly$Comp_Events2 *100 / sum(Monthly$Comp_Events2) ~ c(1:42),
     ylim = c(0, 5),type = "l",col = scales::alpha("blue2", 0.2),pch = 17,xlab = "",
     ylab = "",lty = 3,xaxt = "n",yaxt = "n", bty = "n")

par(new = T)
plot(Monthly$Comp_Events1 *100 / sum(Monthly$Comp_Events1) ~ c(1:42),
     ylim = c(0, 5),type = "l",col = scales::alpha("blue1", 0.2),pch = 18,xlab = "",
     ylab = "",lty = 4,xaxt = "n",yaxt = "n", bty = "n")
legend(1,5,legend = c("Difficult Cases","Comparison Cases (Main)","Comparison Cases (Alternative I)","Comparison Cases (Alternative II)","All Cases"),
       lty = c(1,4,3,2,5), bty = "n", col = c("darkred","blue3","blue1","blue2","palegreen3"), cex = 0.75)
dev.off()


# Visits duration (Figure A4)
pdf(file.path(pd,"outputs","Figure_A4.pdf"))
plot(density(Events$visit_length[Events$Treat==1],na.rm=T), main = "",
     xlim = c(-2,40), ylim = c(0,0.2), col = "darkred", bty = "n", xaxt = "n",
     xlab = "Visit Duration (Minutes)")
par(new=T)
plot(density(Events$visit_length[Events$Control==1],na.rm=T),main = "", ylab = "",lty = 3,
     xlim = c(-2,40), ylim = c(0,0.2), col = "blue1", bty = "n", xaxt = "n",
     xlab = "", yaxt = "n")
par(new=T)
plot(density(Events$visit_length[Events$Control==2],na.rm=T),main = "", ylab = "",lty = 1,
     xlim = c(-2,40), ylim = c(0,0.2), col = "blue2", bty = "n", xaxt = "n",
     xlab = "", yaxt = "n")
par(new=T)
plot(density(Events$visit_length[Events$Control==3],na.rm=T),main = "", ylab = "",lty = 2,
     xlim = c(-2,40), ylim = c(0,0.2), col = "blue3", bty = "n", xaxt = "n",
     xlab = "", yaxt = "n")
axis(1,seq(-5,50,5),pos=0)
par(new=T)
plot(density(Vis$visit_length,na.rm=T),main = "", ylab = "",lty = 5,
     xlim = c(-2,40), ylim = c(0,0.2), col = "palegreen3", bty = "n", xaxt = "n",
     xlab = "", yaxt = "n")
legend(20,0.2,legend = c("Difficult Cases","Comparison Cases (Main)","Comparison Cases (Alternative I)","Comparison Cases (Alternative II)","All Cases"),
       lty = c(1,4,3,2,5), bty = "n", col = c("darkred","blue3","blue1","blue2","palegreen3"), cex = 0.75)
dev.off()

###########################################################################

## Balancing tables (Table 1 and A3)
###########################################################################

# Exclude the index case
Vis <- Vis[Vis$Visits_From_Event != 0, ]

# Sample of visits before the index case (main comparison)
Pre <- Vis[Vis$Visits_From_Event < 0 & Vis$Control %in% c(0,3), ]

# Patient's age and gender
Pre$age <- Pre$year - Pre$birth_date
Pre$Male <- Pre$gender

# Indicator for cancer-related screening test
Pre$Cancer_Test <- (Pre$Mammo + Pre$PSA +Pre$Gastro)

# Balance table (main)
vars <-c("age","Male","ses_ags","TIA","Diabetic","CVD","Hashmana","Cancer_Old",
         "visit_length","Diagnostic_Test","Cancer_Test","ER","Consult","doc_age","doc_gender","doc_exper")

Balance <- NULL
for (v in vars) {
  Pre[[v]] <- Pre[[v]] * 100
  
  if(v == "visit_length"){
    reg <-plm(paste(v, "~Treat", sep = ""),data = Pre[Pre$Visits_From_Event<(-1),],index=c("Real_Event_ID"),model="within")
  }else{
    reg <-plm(paste(v, "~Treat", sep = ""),data = Pre,index=c("Real_Event_ID"),model="within")
  }
  
  co <- coeftest_cluster(reg)
  
  Balance <- rbind (Balance, as.data.frame(
    list(
      "Treat" = (mean(fixef(reg))+co[1,1]),
      "Comparison" = mean(fixef(reg)),
      "Diff" = co[1, 1] ,
      "P_val" = co[1, 4]
    )
  ))
}
rownames(Balance) <- vars
Balance["age", 1:3] <- Balance["age", 1:3] / 100
Balance["ses_ags", 1:3] <- Balance["ses_ags", 1:3] / 100
Balance["doc_age", 1:3] <- Balance["doc_age", 1:3] / 100
Balance["doc_exper", 1:3] <- Balance["doc_exper", 1:3] / 100
Balance["visit_length", 1:3] <- Balance["visit_length", 1:3] / 100
Balance <- round(Balance, 2)

# Export
a <- print(xtable(Balance))
write.table(a,file.path(pd,"outputs","Table_1.txt"))
BalanceM <- Balance

# summarise total numbers
length(unique(Pre$doctor_id))
length(unique(Pre$Real_Event_ID))
length(unique(Pre$Event_ID))-
  length(unique(Pre$Real_Event_ID))
sum(Pre$Treat==1)
sum(Pre$Treat==0)


# Balance table (alternative I comparison)
Pre <- Vis[Vis$Visits_From_Event < 0 & Vis$Control %in% c(0,1), ]
Pre$age <- Pre$year - Pre$birth_date
Pre$Male <- Pre$gender
Pre$Cancer_Test <- (Pre$Mammo + Pre$PSA +Pre$Gastro)
vars <-c("age","Male","ses_ags","TIA","Diabetic","CVD","Hashmana","Cancer_Old",
         "visit_length","Diagnostic_Test","Cancer_Test","ER","Consult","doc_age","doc_gender","doc_exper")

Balance <- NULL
for (v in vars) {
  Pre[[v]] <- Pre[[v]] * 100
  
  if(v == "visit_length"){
    reg <-plm(paste(v, "~Treat", sep = ""),data = Pre[Pre$Visits_From_Event<(-1),],index=c("Real_Event_ID"),model="within")
  }else{
    reg <-plm(paste(v, "~Treat", sep = ""),data = Pre,index=c("Real_Event_ID"),model="within")
  }
  
  co <- coeftest_cluster(reg)
  
  Balance <- rbind (Balance, as.data.frame(
    list(
      "Treat" = (mean(fixef(reg))+co[1,1]),
      "Comparison" = mean(fixef(reg)),
      "Diff" = co[1, 1] ,
      "P_val" = co[1, 4]
    )
  ))
}
rownames(Balance) <- vars
Balance["age", 1:3] <- Balance["age", 1:3] / 100
Balance["ses_ags", 1:3] <- Balance["ses_ags", 1:3] / 100
Balance["visit_length", 1:3] <- Balance["visit_length", 1:3] / 100
Balance["doc_age", 1:3] <- Balance["doc_age", 1:3] / 100
Balance["doc_exper", 1:3] <- Balance["doc_exper", 1:3] / 100
Balance <- round(Balance, 2)
BalanceA1 <- Balance
length(unique(Pre$doctor_id))
length(unique(Pre$Real_Event_ID))
length(unique(Pre$Event_ID))-
  length(unique(Pre$Real_Event_ID))
sum(Pre$Treat==1)
sum(Pre$Treat==0)


# Balance table (alternative II comparison)
Pre <- Vis[Vis$Visits_From_Event < 0 & Vis$Control %in% c(0,2), ]
Pre$age <- Pre$year - Pre$birth_date
Pre$Male <- Pre$gender
Pre$Cancer_Test <- (Pre$Mammo + Pre$PSA +Pre$Gastro)
vars <-c("age","Male","ses_ags","TIA","Diabetic","CVD","Hashmana","Cancer_Old",
         "visit_length","Diagnostic_Test","Cancer_Test","ER","Consult","doc_age","doc_gender","doc_exper")

Balance <- NULL
for (v in vars) {
  Pre[[v]] <- Pre[[v]] * 100
  
  if(v == "visit_length"){
    reg <-plm(paste(v, "~Treat", sep = ""),data = Pre[Pre$Visits_From_Event<(-1),],index=c("Real_Event_ID"),model="within")
  }else{
    reg <-plm(paste(v, "~Treat", sep = ""),data = Pre,index=c("Real_Event_ID"),model="within")
  }
  
  co <- coeftest_cluster(reg)
  
  Balance <- rbind (Balance, as.data.frame(
    list(
      "Treat" = (mean(fixef(reg))+co[1,1]),
      "Comparison" = mean(fixef(reg)),
      "Diff" = co[1, 1] ,
      "P_val" = co[1, 4]
    )
  ))
}
rownames(Balance) <- vars
Balance["age", 1:3] <- Balance["age", 1:3] / 100
Balance["ses_ags", 1:3] <- Balance["ses_ags", 1:3] / 100
Balance["doc_age", 1:3] <- Balance["doc_age", 1:3] / 100
Balance["doc_exper", 1:3] <- Balance["doc_exper", 1:3] / 100
Balance["visit_length", 1:3] <- Balance["visit_length", 1:3] / 100
Balance <- round(Balance, 2)

# Export Table A3
Bal <- cbind(BalanceM[,1:2],BalanceA1[,2],Balance[,2])
a <- print(xtable(Bal))
write.table(a,file.path(pd,"outputs","Table_A3.txt"))

# summarise total numbers
length(unique(Pre$doctor_id))
length(unique(Pre$Real_Event_ID))
length(unique(Pre$Event_ID))-
  length(unique(Pre$Real_Event_ID))
sum(Pre$Treat==1)
sum(Pre$Treat==0)

###########################################################################

## Fix Columns
###########################################################################

Vis$age_sq <- Vis$age^2

Vis$Diagnostic_Test <- Vis$Diagnostic_Test * 100

Vis$Diagnostic_3 <- ((Vis$Blood + Vis$Alt_GPT + Vis$Phosphatase +
                          Vis$Rentgen + Vis$Ultrasound + Vis$Machon)>0) * 100

Vis$Diagnostic_7 <- ((Vis$Blood + Vis$Alt_GPT + Vis$Phosphatase + Vis$Cholesterol +Vis$Triglycerides+Vis$Urine+Vis$TSH+
                        Vis$Ultrasound + Vis$Machon +Vis$Dimut +Vis$Heart_Machon +Vis$Rentgen +Vis$Mammo +Vis$EKG)>0) * 100

Vis$Diagnostic_Number <- ((Vis$Blood + Vis$Alt_GPT + Vis$Phosphatase + Vis$Cholesterol +Vis$Triglycerides+
                             Vis$Ultrasound + Vis$Machon +Vis$Dimut +Vis$Heart_Machon +Vis$Rentgen)) * 100

Vis$Screening <- Vis$Screening_Test * 100

Vis$Lab <- Vis$Lab_Test * 100

Vis$Cancer_Old <- Vis$Cancer_Old * 100

Vis$gender <- Vis$gender * 100

Vis$Cancer_Test <- 100 * (Vis$Mammo + Vis$PSA +Vis$Gastro)

Vis$Hospitalizations_Last_6_Months <- Vis$Hospitalizations_Last_6_Months * 100

Vis$Visits_Last_6_Months <- Vis$Visits_Last_6_Months * 100

Vis$Consult <- (Vis$Consult>0) * 100

Vis$ER <- (Vis$ER) * 100

Vis$Drug <- (Vis$Drug) * 100

###########################################################################

## Analysis
###########################################################################

# Figure 1 - Dynamic DD
dynamic_dd(outcome = "Diagnostic_Test",
                  Label =  "Rate of Referrals to Tests (p.p.)",
                  min1 =  -4,
                  max1 =  4)

# Figure A5 and 2 - PS calibration plot and effects along the estimated PS
PS_DDD(outcome="Predicted_xg_fit",CON = 3,v1 = 8)

# Table A1 
cancers <- summarise(group_by(Events[Events$Treat==1,],Cancer),total=n(),Mortality=mean(Death_4_Years,na.rm=T)) # Fix to 1
stargazer(cancers,type="latex",out="Emotions/outputs/Table_A1.txt",summary = F)

# Table A2- referrals frequencies
V <- c("Blood","Alt_GPT","Phosphatase","Cholesterol","Triglycerides",
       "Urine","TSH","AST_GOT","VitaminB12","Ferritin",
       "Rentgen","Ultrasound","Machon","Dimut","Heart_Machon",
       "Mammo","EKG","Density","Dufler","Heart_Eco")
for(v in V){
  print(mean(Vis[[v]]))
}

# Table A4 - placebo
DD_main(outcomes = c("age","Cancer_Old","gender","Visits_Last_6_Months","Hospitalizations_Last_6_Months"),
        names = c("Patient Age (Years)","Previous Cancer Diagnosis (Percent)",
                  "Male Patient (Percent)",
                  "Any Physician Visit in 6 Months (Percent)",
                  "Any Hospital Admission in 6 Months (Percent)"),
        CON = 3, name = "Table_A4",
        v1=8)

# Table 2 - DD, main outcomes
DD_main(outcomes = c("Diagnostic_Test","Drug","Consult","ER","Cancer_Test","visit_length"),
         names = c("Test Referral","Drug Prescription","Specialist Referral","ER Referral","Cancer Screening Test Referral","Visit Duration (minutes)"),
         CON = 3, 
         v1 = 8,
         name = "Table_2")

# Table 3 - DDD, testing
Vis$Above_56_Doc <- as.numeric(Vis$doc_age>56)
Vis$Deadly <- Vis$Cancer %in% c("Pancreas","Liver","Stomach","Lung","Leukemia","Myeloma","Kidney",
                                "Bladder","Colon","Lymphoma","Rectum","Ovary")
DDD(outcome = "Diagnostic_Test",
    Z = c("Above_56_Doc","doc_gender","Exposure_5","Above_62_Event","Deadly","Death_4_Years"),
    names = c("Age > 56","Male","Exposure >5","Age > 62","High-Risk Cancer","Died Within 4 Years"),
    v1 = 8,
    CON = 3,
    name = "Table_3")


# Table A5 - DD with alternative comparison groups
DD_main_alternative(outcomes = c("Diagnostic_Test","Drug","Consult","ER","Cancer_Test","visit_length"),
             v1 = 8,
             names = c("Test Referral","Drug Prescription","Specialist Referral","ER Referral","Cancer Screening Test Referral","Visit Duration (minutes)"),
            name = "Table_A5")

# Table A6 - DD with controls
DD_main_controls(outcomes = c("Diagnostic_Test","Drug","Consult","ER","Cancer_Test","visit_length"),
                  v1 = 8,
                  names = c("Test Referral","Drug Prescription","Specialist Referral","ER Referral","Cancer Screening Test Referral","Visit Duration (minutes)"),
                  CON = 3,
                 name = "Table_A6")

# Table A7 - DD with Alternative outcomes
DD_main(outcomes = c("Lab","Screening","Diagnostic_7","Diagnostic_3","Diagnostic_Number"),
         names = c("Top 5 Lab Tests","Top 5 Imaging and Other Tests","Top 7 Tests","Top 3 Tests","Number of Tests"),
         CON = 3, 
         v1 = 8,
        name = "Table_A7")

# Start/Finish Day Timing Analysis (Table A9)
days <- distinct(Vis,Event_ID,.keep_all = T)
days3 <- days[days$Control %in% c(0,1), ]
days3$Finish_The_Day_Time <- floor(days3$Finish_The_Day_Time)
days3$Start_The_Day_Time <- floor(days3$Start_The_Day_Time)
lmstart3 <- plm(Start_The_Day_Time ~ Treat, data = days3,index = c("Real_Event_ID"),model = "within")
lmend3 <- plm(Finish_The_Day_Time ~ Treat, data = days3,index = c("Real_Event_ID"),model = "within")
stargazer(lmstart3,lmend3,type="latex",out="Emotions/outputs/Table_A8.txt",
          add.lines = list(c("Baseline",round(mean(days3$Start_The_Day_Time[days3$Treat==0]),2),
                             round(mean(days3$Finish_The_Day_Time[days3$Treat==0]),2)),
                           c("N",round(length(lmstart3$residuals),
                                       length(lmend3$residuals)))))

###########################################################################
